%% Header
% in unit square, made specifically for a square number of plants
% and variation in plant location amenities
% can set various parameter values in following section 
clc;
clear;
addpath Functions;


tic

Qs = [3,4,5,7,9,10,12,13,15,17,18,20,22,23,25];
l = length(Qs);
for idx=1:l
%% params
filename = 'N144_q05_deltapoint45.mat';
delta=9;
q=Qs(idx);

epsilon=2;
theta=2;
xi=delta^2;
grid_sz=6; % # of points (s) to make in square, grid_sz+1=dimension

Ds=zeros(grid_sz+1,grid_sz+1); % set up demand at s locations
for i=0:grid_sz
    for j=0:grid_sz
        x=i/grid_sz;
        y=j/grid_sz;
        Ds(i+1,j+1)=-1.5*sin((3*pi)*x*(pi*y))+1.5;
    end
end

sigma=0.45; %low will encourage more plants

side=6; %number of possible plant locations along 1 side
plant_num=side^2; % number of possible plant locations total
plant_locs=zeros(plant_num,2);
for i=1:side
    for j=1:side
        plant_locs(j+(i-1)*side,:)=[1/(side*2)+(i-1)/side, 1/(side*2)+(j-1)/side ];
    end
end

BoW=ones(length(plant_locs),1)*1; % uniform BoW

R=zeros(length(plant_locs),1); % set up rents
for k=1:length(plant_locs)
    x=plant_locs(k,1);
    y=plant_locs(k,2);
    R(k)=1.5*sin((3*pi*x)*(pi*y))+1.5;
end




% make struct to pass into AER
parameters=struct('epsilon', epsilon,...
    'theta', theta,...
    'xi', xi,...
    'q',q,...
    'sigma',sigma,...
    'grid_sz',grid_sz,...
    'plant_locs',plant_locs,...
    'D',Ds,...
    'BoW',BoW,...
    'R',R,...
    'delta',delta);

%% set up vector to make decisions and run
unsure=0.5*ones(length(plant_locs),1);

[choice, profits]=AER(unsure, 'calc_pi', true, parameters, true);
disp(strcat('Number of Operational Plants:',num2str(sum(choice))));
tot_rev=profits-sum(R.*choice);
N_init=sum(choice);

end
toc


